#!/usr/bin/perl -w
#-----------------------------------------------------------+
#                                                           |
# RANDOMIZE SEQ READS                                       |
#                                                           |
#-----------------------------------------------------------+
#                                                           |
#  AUTHOR: James C. Estill                                  |
# CONTACT: JamesEstill @ gmail.com                          |
# STARTED: 08/23/2006                                       |
# UPDATED: 08/25/2006                                       |
# DESCRIPTION:                                              | 
#  For a directory containing subdirs with the zipped data  |
#  downloaded from NCBI, this will unzip, untar and move    |
#  all *.fasta, *.scf, *.qscore, *.reversecom  files to     |
#  the parent dir. A number of sequences will be randomly   |
#  selected from the dir, and the *.fasta and *.qscore      |
#  files sent to concantenated files.                       |
#                                                           |
#-----------------------------------------------------------+

=head1 INCLUDED MODULES
External modules used by the program.
=cut
#----------------------------+
# INCLUDED MODULES           |
#----------------------------+
use DBI;                      # Allows connection to MySQL database
use Bio::SeqIO;
use Bio::Seq;
use Bio::Tools::Lucy;         # Use some Lucy 

=head1 VARIABLES
Program variables
=cut     
#-----------------------------+
# USER SET VARIABLES          |
#-----------------------------+
my $BaseDir = "/home/jestill/projects/maize/RAND/"; 
my $SeqTable = "tblRandMaize";
my $GetScf = "true";          # Include the scf files in the outputdir

#----------------------------+
# VARIABLES                  |
#----------------------------+
my $ProcNum = "0";
my $RecNum = "0";             # The absolute record number
my $DbUserName = "jestill";   # Username for MySQL Database
my $DbName = "dbRandMaize";   # DB Name for MySQL Database
my $DbUserPassword;           # Password for MySQL Database    

#-----------------------------+
# GET THE USER PASSWORD FROM  |
# THE COMMAND LINE            |
#-----------------------------+

print "\nPassword for $DbUserName\n";
system('stty', '-echo') == 0 or die "can't turn off echo: $?";
$DbUserPassword = <STDIN>;
system('stty', 'echo') == 0 or die "can't turn on echo: $?";
chomp $DbUserPassword;

=head1 DATABASE CONNECT & FILE IO
Connect to the database
=cut
#-----------------------------+
# DATABASE CONNECTION         |
#-----------------------------+
my $dbh = DBI->connect("DBI:mysql:database=$DbName;host=localhost",
		       $DbUserName, $DbUserPassword,
		       {'RaiseError' => 1});

#-----------------------------+
# CREATE THE SEQ DATA TABLE   |
#-----------------------------+
&CreateDBTable ( $SeqTable );

=head1 MAIN BODY
Main body of the program.
=cut
#-----------------------------------------------------------+
# MAIN BODY                                                 |
#-----------------------------------------------------------+

# BASE DIR THAT ALL OUTPUT WILL BE IN
my @BacDirs = ("AC145389", # FASTA FILES MISSING FROM NCBI DOWNLOAD
	       "AC146795",
	       "AC147505", 
	       "AC148082",
	       "AC148160",
	       "AC145726",  
	       "AC146813",  
	       "AC147518", 
	       "AC148093", 
	       "AC148163",
	       "AC145729",  
	       "AC147015",
	       "AC147602",  
	       "AC148096",
	       "AC148169",
	       "AC146764",  
	       "AC147504",
	       "AC147776",
	       "AC148112",
	       "AC148394",);

foreach $IndBacDir (@BacDirs)
{
    print "PROCESSING ".$IndBacDir."\n";
    $ProcNum++;

    #if ($ProcNum == 2){exit;}     # Only process first rec for test
    #-----------------------------+
    # EXTRACT SEQUENCE DATA FROM  |
    # THE GENBANK TRACE DOWNLOAD  |
    #-----------------------------+
    
    &UnzipDir ($BaseDir.$IndBacDir."/MIPS_trace.tar.gz",
	       $BaseDir.$IndBacDir."_Out");
    
    &RandFasta ( $BaseDir.$IndBacDir."_Out/",
		 "96",
		 $BaseDir.$IndBacDir."_Rand",
		 $IndBacDir );

    
    my $ClearDirCmd = "rm -rf ".$BaseDir.$IndBacDir."_Out/";
    system ($ClearDirCmd);
    
}

#-----------------------------+
# PROCESS HAS COMPLETED       |
#-----------------------------+
print "\a";
print "THE PROCESS IS COMPLETE\n";

exit;

sub UnzipDir
{
#-----------------------------+
# UNTAR ENTIRE DIR SUBTREE TO |
# A SINGLE OUTPUT DIR         |
#-----------------------------+
# The dirs from the original directory tree will
# still be present
    
    my $InFile = $_[0];        # The gzipped and tar file to unzip
    my $OutDir = $_[1];
    #my $WorkDir = $_[2] || $_[1];
    
    # DEBUG/CHECK VARS
    print "IN FILE:\n".$InFile."\n";
    print "OUT DIR:\n".$OutDir."\n";
    print "WORD DIR:\n".$OutDir."\n";
    
    # CREATE THE OUTPUT DIR IF NEEDED
    mkdir $OutDir, 0777 unless (-e $OutDir);

    # Only try to do the work if the file exists
    if (-e $InFile)
    {
	#-----------------------------+
	# UNTAR THE FILE              |
	#-----------------------------+
	print "UNTAR the file $InFile\n";
	my $TarCmd = "tar -C ".$OutDir." -zxf ".$InFile;
	system ($TarCmd);

	#-----------------------------+
	# MAKE THE OUT DIR THE CURRENT|
	# WORKING DIR                 |
	#-----------------------------+
	chdir ($OutDir) ||
	    die "Can not change to directory\n".$OutDir."\n";

	#-----------------------------+
	# USE THE FIND COMMAND TO MOVE|
	# THE FILES                   |
	#-----------------------------+
	print "MOVING THE FILES TO PARENT\n";
	#my $MoveCmd = "find . -type f -exec mv -i \'{}\' . \\\;";
	# The f flag will force overwrite of any existing files
	my $MoveCmd = "find . -type f -exec mv -if \'{}\' . \\\;";
	print $MoveCmd."\n";
        system ($MoveCmd);

    }else{
	# SOUND THE ALARM AND TELL THE USER THE FILE
	# DOES NOT EXIST
	print "\a"; 
	print "INFILE DOES NOT EXIST\n".$InFile."\n";
	exit;
    }
}


sub RandFasta
{
#-----------------------------+
# RANDOMLY SELECT SEUQNECE    |
# FILES                       |
#-----------------------------+
# Uses a search for the fasta file to build a base name
# that other files will be searched for

    my $SearchDir = $_[0];     # The directory to search
    my $numSel = $_[1];        # The number of reads to select
    my $BaseOut = $_[2];       # The base path for the output files
                               # this should include the dir and
                               # the base name for building out files.
    my $BacID = $_[3];         # The unique id for the BAC
    my $i;
    my @SortAry;

    #-----------------------------+
    # OPEN OUT FILES AND DIR      |
    #-----------------------------+
    opendir(DIR, $SearchDir);
    my @FastaFiles = grep(/\.fasta$/,readdir(DIR));
    closedir(DIR);
    my $numFiles = @FastaFiles;
    my $numReps = $numFiles - 1;

    # CONCATENATED FASTA FILE
    my $FastCatFile = $BaseOut.".fasta";
    open (FASTOUT, ">".$FastCatFile) ||
	die "Can not open:\n".$FastCatFile."\n";
    
    # CONCATENATED QSCORE FILE
    my $QScoreCatFile = $BaseOut.".qscore";
    open (QUALOUT, ">".$QScoreCatFile) ||
	die "Can not open:\n".$QScoreCatFile."\n";
    
    #-----------------------------+
    # THROW AN  ERROR IF TOO MANY |
    # FILES ARE SELECTED          |
    #-----------------------------+
    if ($numSel > $numFiles)
    {
	print "ERROR!!\n";
	print "The number of files to select: ".$numSel."\n";
	print "Exceeds the number of files in the directory: ".
	    $numFiles."\n";
	exit;
    }

    #-----------------------------+
    # LOAD DATA TO A 2D SORTABLE  |
    # ARRAY                       |
    #-----------------------------+
    for ( $i=0; $i <= $numReps; $i++)
    {
	my $FileName = $FastaFiles[$i] ||
	    "UNDEF";
	
	$SortAry[$i][0] = $FileName;
	$SortAry[$i][1] = rand();
    }

    #-----------------------------+
    # SORT THE FASTAFILES BY THE  |
    # RANDOM NUMBER               |
    #-----------------------------+
    @SortedAry = sort {$a->[1] <=> $b->[1]} @SortAry;


    #-----------------------------+
    # FOR THE TOP NUMSELECTED RECS|
    # IN THE SORTED ARRAY         |
    #-----------------------------+
    for ( $i=0; $i < $numSel; $i++)
    {
	$RecNum++;                
	#-----------------------------+
	# GET VARIABLES               |
	#-----------------------------+
	# Fetching the vars to see how it worked
	my $FileName = $SortedAry[$i][0] ||
	    "ERROR";
	my $RandNum = $SortedAry[$i][1] ||
	    "ERROR";

	#-----------------------------+
	# SET VARIABLE SCOPE          |
	#-----------------------------+
	my $BaseInName;                # The base name to build infile names

	#-----------------------------+
	# GET BASE NAME               |
	#-----------------------------+
	if ($FileName =~ m/(.*)\.fasta/)
	{
	    $BaseInName = $1;
	}else{
	    $BaseInName = "MATCH ERROR";
	    print "Error matching the base fasta input name\n".
		$FileName."\n";
	    exit;
	};
		
	#-----------------------------+
	# COPY THE INPUT FASTA FILE TO| 
	# THE CONCATENATED FILE       |
	#-----------------------------+
	my $FastFilePath = $SearchDir.$BaseInName.".fasta";
	open (FASTIN, "<".$FastFilePath) ||
	    die "Can not open input file\n".$FastFilePath."\n";
	while (<FASTIN>)
	{
	    #-----------------------------+
	    # IF THIS IS THE HEADER LINE  |
	    # THEN APPEND A UNIQUE ID NUM |
	    # AS A PREFIX TO THE NAME     |
	    #-----------------------------+
	    my $FirstChar = substr ($_, 0, 1);
	    if ($FirstChar =~ ">")
	    {
		my $TrimRec = substr ($_, 1);
		print FASTOUT ">".$RecNum."|".$TrimRec;
		#print FASTOUT $RecNum."|".$_;
	    }else{
		print FASTOUT $_;
	    }
	}
	close FASTIN;
	

	#-----------------------------+
	# LOAD SEQUENCE DATA TO THE   |
	# DATABASE                    |
	#-----------------------------+
	&UploadSeq ( $BacID,
		     $FastFilePath,
		     $RecNum);
	
	#-----------------------------+
	# COPY THE QUALITY SCORE  FILE|
	# TO THE CONCATENATED FILE    |
	#-----------------------------+
	my $QScorePath = $SearchDir.$BaseInName.".qscore";
	open (QUALIN, "<".$QScorePath) ||
	    die "ERROR!! Can not open input file path:\n".$QScorePath."\n";
	while (<QUALIN>)
	{
	    my $FirstChar = substr ($_, 0, 1);
	    if ($FirstChar =~ ">")
	    {
		my $TrimRec = substr ($_, 1);
		print QUALOUT ">".$RecNum."|".$TrimRec;
	    }else{
		print QUALOUT $_;
	    }
	    
	    #print QUALOUT $_; 
	}
	close QUALIN;


    } # End of for each file selected $i is 1 to numSelected
    
    print "\n\n".$numSel." files successfully selected from\n".
	$numFiles." in the dir.\n";

    close FASTOUT;
    close QUALOUT;
    
    #-----------------------------+
    # CLEAN THE FASTA OUTPUT USING|
    # THE TIGR PROGRAM LUCY       |
    #-----------------------------+
    my $LucyFasta = $BaseOut.".lucy.fasta";
    my $LucyQual = $BaseOut.".lucy.qscore";
    &SeqClean ($FastCatFile, $QScoreCatFile, 
	       $LucyFasta, $LucyQual);
 

    #-----------------------------+
    # UPDATE DATABASE WITH INFO   |
    # FROM LUCY CLEANED FLIE      |
    #-----------------------------+
    &LucyUpdateDb ($LucyFasta,
		   $LucyQual);
   
}

sub SeqClean 
{
#-----------------------------+
# CLEAN SEQUENCE FILES USING  |
# THE LUCY PROGRAM            |
#-----------------------------+
# LUCY available from TIGR
# http://www.tigr.org/software/index.shtml

    my $LucFastaIn = $_[0];    # FASTA file to send to LUCY
    my $LucQualIn = $_[1];     # Quality file to send to LUCY
    my $LucFastaOut = $_[2];   # Lucy tagged FASTA out
    my $LucQualOut = $_[3];    # Lucy tagged quality output
 
    my $LucVecFile = "/home/jestill/Apps/lucy/lucy-1.19p/PUC19";
    my $LucSpliceFile = "/home/jestill/Apps/lucy/lucy-1.19p/PUC19splice";
   
    #-----------------------------+
    # RUN LUCY TO CLEAN THE       |
    # FASTA FILES                 |
    #-----------------------------+
    # Would need to add an info file
    my $LucyCmd = "lucy -v ".$LucVecFile." ".$LucSpliceFile." ".
	"-output ".$LucFastaOut." ".$LucQualOut." ".
	$LucFastaIn." ".$LucQualIn;
    
    print $LucyCmd."\n";           # Show the command
    system ($LucyCmd);
}

sub CreateDBTable
{
    my $CreateTableName = $_[0];   
 
    # DROP THE TABLE it already exists
    if (&does_table_exist( $dbh, $CreateTableName ))
    {
	$dbh->do("DROP TABLE ".$CreateTableName);
    }

    $qryCreateTable = "CREATE TABLE ".$CreateTableName." (".
	" rownum INT(10) NOT NULL AUTO_INCREMENT,".
	" bac CHAR(10),".     # Bac sublcone was derived from
	" src CHAR(4),".      # Source of subclone (UGA/TAMU)
	" seq_run char(5),".  # Sequence run (# Of Seq Run)
	" seq_set char(10),". # Sequence set (Random/Targeted)
	" plt_addr char(5),". # Sequencing Plate Address
	" primer char(10),".  # Primer used for seqeuncing
	" species char(20),". # Species the seq is derived from
	" id char(40),".      # Long unique ID for the Seq
	" seq MEDIUMTEXT,".   # Sequence for the read
	" spec_num char(8),". # Number within species, used for XY plot
	" seq_len char(5),".  # Length of the seq
	" abt MEDIUMTEXT,".   # Description of the sequence
	" luc_start INT(5),". # Start of LUCY good qual sequence
	" luc_end INT(5),".   # End of the LUCY good qual sequence
	" luc_len INT(5),".   # Length of the LUCY good qual sequence
	" luc_seq MEDIUMTEXT,". # LUCY CLEAN SEQUENCE
	" KEY (rownum))";

    $dbh->do($qryCreateTable);

}

sub does_table_exist
{
    my ($dbh,$whichtable) = @_;
    my ($table,@alltables,$found);
    @alltables = $dbh->tables();
    $found = 0;
    foreach $table (@alltables) {
	$found=1 if ($table eq "`".$whichtable."`");
    }
    return $found;
}


sub UploadSeq


{
#-----------------------------+
# LOAD SEQUENCE DATA FROM THE |
# HIGH QUALITY READ TO THE    |
# DATABASE                    |
#-----------------------------+
    print "UPLOADING DATA TO DB\n";
    my $BAC = $_[0];
    my $SeqFile = $_[1];
    my $UniqPre = $_[2];
#    my $SeqTable = $_[2];
    my $SeqFormat = "fasta";
    my $Src = "NCBI";          # This data is from NCBI tracefile
    my $SeqRun = "NA";
    my $Species = "maize";
    my $SpecNum = "1";
    my $inseq = Bio::SeqIO->new(-file   => "<$SeqFile",
				-format => $SeqFormat ) 
	|| die "ERROR!! Can not open:\n $SeqFile\n";    


    while (my $seq = $inseq->next_seq) 
    {
	$SeqUniqueId = $seq->primary_id;
	# Add prefix unique id to the name
	$SeqUniqueId = $UniqPre."|".$SeqUniqueId;

	#------------------------------+
	# UPLOAD INFO TO THE DATABASE  |
	#------------------------------+
	$UploadData = "INSERT IGNORE INTO ".$SeqTable.
	    " ( bac, src, seq_run,".
	    " species, id, spec_num, seq, seq_len, abt)".
	    " VALUES (".
	    " '".$BAC."',".
	    " '".$Src."',".
	    " '".$SeqRun."',".
	    " '".$Species."',".
	    " '".$SeqUniqueId."',".
	    " '".$SpecNum."',".
	    " '".$seq->seq."',".
	    " '".$seq->length."',".
	    " '".$seq->desc."'".
	    " )";
	#print "SQL:".$UploadData."\n"; # Only show this for Error Checking
	
	$dbh->do($UploadData);
    }
}

sub LucyUpdateDb 
{
#-----------------------------+
# UPDATE SEQUENCE DATA WITH   |
# CLEAN SEQ DATA FROM LUCY    |
#-----------------------------+

# NOTE QUALITY APPEARS TO BE WRAPPING THE
# RETURN CHARACTER, PROBABLY WRITTEN 
# FOR WINDOWS, NOT UNIX LINE ENDINGS

    print "ADDING LUCY DATA TO THE DB\n";
    
    my $LucyFastaFile = $_[0];
    my $LucyQualFile = $_[1];

    my @LucyParams = ('seqfile' => $LucyFastaFile, 
		      #'qualfile' => $LucyQualFile,
		      'lucy_verbose' => 1);

    my $LucyObj = Bio::Tools::Lucy->new(@LucyParams);
    
    $LucyNames = $LucyObj->get_sequence_names();
    
    foreach $IndName (@$LucyNames) 
    {
	#-----------------------------+
	# WORK THROUGH PARSING LUCY   |
	#-----------------------------+
	my $ClearLen = $LucyObj->length_clear($IndName);
	my $ClearStart = $LucyObj->start_clear($IndName);
	my $ClearEnd = $LucyObj->end_clear($IndName);
	my $LucSeq = $LucyObj->sequence($IndName);
	my $LucSeqId = $IndName;
	
	print $LucSeqId."\n";
	print "\tSTART\t".$ClearStart."\n";
	print "\tEND\t".$ClearEnd."\n";
	print "\tLEN\t".$ClearLen."\n";

	my $qryAddLucy = "UPDATE ".$SeqTable.
	    " SET ".
	    " luc_start=\'".$ClearStart."\',".
	    " luc_end=\'".$ClearEnd."\',".
	    " luc_len=\'".$ClearLen."\',".
	    " luc_seq=\'".$LucSeq."\'".
	    " WHERE id=\'".$LucSeqId."'";

	$dbh->do($qryAddLucy);

    }

    #-----------------------------+
    # IF A STANDARD ERROR FILE WAS|
    # PRODUCED YOU CAN FETCH WHY  |
    # FILES WERE REJECTED         |
    #-----------------------------+ 


}

=head1 HISTORY
=cut
#-----------------------------------------------------------+
# HISTORY                                                   |
#-----------------------------------------------------------+
# 08/23/2006
# - Program started
# - Added UnzipDir subfunction
# - Added RandFasta subfunction
#
# 08/24/2006
# - Changed RandFasta to select all files, not just the FASTA
#   files for output to a single set of cat files
#
# 08/25/2006
# - Worked more on the RandFasta function to take input from other
#   qual files etc if avaialable. 
# - Added a SeqClean function step using LUCY
# - Added CreateDBTable subfunction
# - Added does_table_exist subfunction
# - Added UploadSeq subfunction
# 
# 08/28/2006
# - Finishing the UploadSeq subfunction, mostly copied code 
#   from ParseASGR4.pl
# - Added LucyUpdateDb subfunction
#   NOTE: Tried to parse quality information from LUCY output
#         but BioPERL parser not working. Appears to be ingnoring
#         line endings when parsing quality info.
#
# 08/29/2006
# - Modified the RandSeq output to include a unique incremental
#   numerical ID. This required changing output to the 
#   concatenated FASTA and QUALITY files as well as the function
#   to upload data to the db.
# 
#-----------------------------------------------------------+
# TO DO                                                     |
#-----------------------------------------------------------+


